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INTRODUCTION 


The  computational  cost  associated  with  the  analysis,  design,  and  optimization  of  the  transient 
behavior  of  a  nonlinear  mechanical  system  can  be  very  high.  Such  problems  are  typically  solved 
with  time  marching  schemes,  via  a  temporal  finite  difference  of  the  set  of  ordinary  differential 
equations  (which  are  commonly  formulated  with  spatial  finite  element  methods).  Implicit  time 
marching  schemes  are  preferred  for  long-duration  structural  dynamics  problems1,  which  require 
the  solution  of  a  system  of  equations  many  times  within  each  time  step  for  a  nonlinear  system  (to 
drive  the  residual  to  zero  with,  for  example,  a  Newton-Raphson  update  loop').  Conservation  of 
energy  within  each  time  step  of  an  implicit  scheme  can  be  difficult  to  enforce  under  the  presence 
of  strong  nonlinearities  ,  potentially  leading  to  instabilities.  Periodically-actuated  structures  may 
take  many  cycles  for  the  response  to  set  up  into  a  time-periodic  state;  time  marching  schemes  can 
be  inefficient  if  the  transients  decay  slowly4. 

Finally,  in  the  event  that  a  large  number  of  design  parameters  are  of  interest  for  system 
optimization,  analytical  sensitivities  of  a  metric  of  the  unsteady  nonlinear  response  with  respect 
to  these  variables  must  be  provided  to  the  optimization  for  a  gradient-based  search.  It  is  well- 
known  that  for  design  studies  with  a  larger  number  of  design  variables  than  constraints,  adjoint 
methods  are  less  expensive  than  direct  methods5.  The  latter  directly  computes  the  derivative  of 
the  system  response  with  respect  to  the  design  variables,  while  the  former  skips  this  (typically) 
superfluous  quantity  in  favor  of  an  adjoint  vector,  which  is  independent  of  the  design  variables. 
The  adjoint  method,  however,  is  notoriously  cumbersome  for  nonlinear  time  marching 
problems6.  The  differential  equation  for  the  adjoint  vector  is  provided  with  terminal  conditions, 
and  must  be  integrated  in  reverse  :  the  system  response  and  the  adjoint  cannot  be  computed 
simultaneously,  and  the  storage  costs  can  be  large. 

One  potential  solution  to  all  of  the  aforementioned  issues  is  a  time -periodic  scheme.  These 
methods  approximate  the  time  derivative  operators  in  the  set  of  ordinary  differential  equations  as 
a  set  of  algebraic  equations,  which  can  then  be  solved  to  compute  the  entire  time  response 
simultaneously.  Furthermore,  a  time-periodic  condition  can  be  enforced  within  the  system  to 
compute  only  the  time-periodic  response,  bypassing  the  initial  transients  completely.  Having 
solved  the  nonlinear  algebraic  system  of  equations  for  the  time -periodic  response,  the  adjoint 
vector  (and  subsequent  design  derivatives)  are  easily  and  inexpensively  computed  '  .  Three 
general  methods  are  commonly  used  for  such  a  time-periodic  analysis:  harmonic  balance 
methods,  cyclic-implicit  methods,  and  spectral  element  methods.  The  harmonic  balance  method 
involves  approximating  each  structural  degree  of  freedom  as  a  Fourier  series,  and  then  solving 
for  the  coefficients  associated  with  each  harmonic10,  or  recasting  the  problem  in  the  time 
domain11,  among  a  variety  of  approaches.  The  cyclic -implicit  method  uses  a  time  marching 
scheme  over  an  actuation  cycle  to  produce  the  relationship  between  the  response  at  adjacent  time 
steps,  and  then  sets  the  system  response  at  the  final  step  identical  to  that  at  the  first  step  . 
Finally,  the  spectral  element  method  discretizes  the  time  domain  into  elements13'  14,  with  a 
polynomial  approximation  of  the  solution  within  each  spectral  element4. 

The  harmonic  balance  method  may  suffer  from  the  global  nature  of  the  Fourier 
approximation  in  terms  of  computational  efficiency15  for  systems  with  a  large  range  of  periodic 
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content.  Furthermore,  excessive  nonlinearity  can  cause  the  higher  harmonics  to  corrupt  the 
lower  harmonics  (aliasing),  producing  non-physical  solutions  and  non-convergence  of  the 
Fourier  series11;  the  harmonic  balance  method  will  not  be  considered  in  this  work.  Of  the  two 
remaining  methods  (cyclic-implicit  and  spectral  elements),  two  substantial  issues  are 
encountered  upon  implementation.  The  first  is  that  the  Jacobian  of  either  method  becomes 
relatively  ill-conditioned  as  the  time -periodic  assumption  is  enforced4,  leading  to  convergence 
issues  within  a  Newton-Raphson  loop  if  the  nonlinearities  are  strong.  The  second  issue  is  that 
the  size  of  the  Jacobian  can  be  extremely  large:  the  product  of  the  number  of  spatial  and 
temporal  degrees  of  freedom.  The  computational  cost  of  solving  this  system  can  be  drastically 
alleviated  by  projecting  the  system  of  ordinary  differential  equations  onto  a  reduced  basis,  prior 
to  conversion  to  a  set  of  algebraic  equations  via  either  spectral  elements  or  the  cyclic-implicit 
method. 

Many  basis  functions  can  be  reasonably  used,  though  POD  modes  (proper  orthogonal 
decomposition)  are  the  most  popular  for  nonlinear  systems,  and  will  be  used  here.  POD  is  a 
linear  method  whose  bases  minimize  the  Euclidean  norm  of  the  distance  between  the  full  order 
response  and  the  reduced  generalized  coordinates;  as  such,  for  a  given  number  of  modes  no  other 
basis  will  be  able  to  capture  as  much  information16.  In  general,  very  few  POD  modes  (relative  to 
the  number  of  spatial  degrees  of  freedom)  are  needed  for  an  accurate  reconstruction  ,  though  the 
modes  must  be  computed  by  obtaining  samples  (snapshots)  of  the  full  order  nonlinear  transient 
analysis.  Data  from  the  time-periodic  methods  ought  not  be  used  for  snapshots,  as  these  methods 
compute  the  entire  response  simultaneously,  rendering  the  model  reduction  obsolete. 
Alternatively,  implicit  integration  can  be  run  through  a  training  period,  the  POD  modes  can  be 
computed  to  construct  the  reduced  order  model  (ROM),  which  can  then  be  converted  into  a  set  of 
algebraic  system  of  equations  for  time-periodic  analysis.  Having  done  this,  the  reduced  cyclic 
adjoint  vector  is  computed  for  design  sensitivities. 

The  work  presents  a  comparative  analysis  of  six  methods  for  computation  of  system  response 
and  adjoint  sensitivities  of  a  nonlinear  structure  under  periodic  actuation:  time  marching  implicit 
integration,  the  cyclic-implicit  method,  and  the  spectral  element  method,  with  each  method 
considered  both  with  and  without  POD-based  model  reduction.  A  brief  derivation  is  provided 
for  each,  for  both  system  response  and  adjoint  sensitivity  computations.  The  salient  issues 
pertaining  to  ROM  construction  are  provided  as  well.  The  system  under  consideration  is  a  two- 
dimensional  elastic  beam  periodically  actuated  at  its  base  through  a  prescribed  rotation.  The 
computational  cost  of  structural  optimization  of  the  beam  is  provided  for  cases  with  weak  and 
strong  geometric  nonlinearities,  with  a  low  (i.e.,  sinusoidal)  and  high  harmonic  content,  for  each 
of  the  six  methods. 
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NONLINEAR  BEAM  DYNAMICS 


This  work  is  concerned  with  the  nonlinear  deformations  of  a  planar  elastic  beam,  periodically 
actuated  at  its  base,  as  seen  in  Figure  1.  Two  different  coordinate  systems  must  be  used  to 
develop  the  equations  of  motion:  a  floating  frame  which  rotates  with  the  beam  (x-y),  and  is 
attached  to  it  base,  and  a  fixed  inertial  frame  (X-Y).  The  beam  is  discretized  into  finite  elements, 
and  each  node  is  given  three  degrees  of  freedom:  axial  displacement,  transverse  displacement, 
and  rotation.  The  beam  is  considered  to  be  clamped  at  its  base,  in  the  floating  frame,  and  the 
deformations  are  computed  in  the  floating  frame  as  well.  The  length  of  the  beam  is  1  m,  the 
square  cross-sectional  area  is  4  cm2,  the  density  is  50  kg/m3,  and  the  elastic  modulus  is  600  MPa. 
The  prescribed  rotation  is: 


6(t)  = 


Tt 

4  •  sin_1(K) 


•  sin 


(1) 


The  amplitude  is  fixed  at  45°,  and  K  lies  between  0  and  1.  At  the  lower  limit,  9  becomes 
sinusoidal,  while  at  the  upper  limit  9  approaches  a  triangular  waveform  .  For  this  work,  K 
values  of  0.01  and  0.99  are  studied  (as  seen  in  Figure  2),  with  actuation  frequencies  (co)  of  15 
and  40  rad/s  (the  latter  of  which  imposes  strong  geometric  nonlinearities  upon  the  beam,  from 
inertial  loads). 

The  equations  of  motion  for  the  beam  are: 


[M]  •  {u}  +  [C]  •  {u}  +  {P}  =  {F} 


[K] 


d{P } 
3{u} 


(2) 


where  the  deformation  vector  { u }  is  as  defined  above  (measured  in  the  body-fixed  coordinate 
system),  [M]  is  a  consistent  mass  matrix,  [C]  is  a  damping  matrix  composed  of  skew-symmetric 
Coriolis  terms  and  symmetric  damping  terms  (Rayleigh  damping,  given  by  10- [M]),  the  internal 
force  { P }  is  a  vector  of  internal  forces  (a  nonlinear  function  of  the  deformation  { u } )  and  { F }  is 
an  external  force  vector  due  to  the  rigid  body  inertial  motions.  Though  not  explicitly  included  in 
Eq.  2,  the  tangent  stiffness  [K]  is  needed  for  a  Newton-Raphson  iteration  scheme.  Inertial 
portions  of  Eq.  2  are  computed  with  standard  multibody  dynamics  techniques,  while  elastic 
portions  are  computed  with  a  corotational  method19.  The  solution  of  this  set  of  nonlinear 
equations  of  motion  will  provide  the  beam  deformation  { u }  as  a  function  of  time. 
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Figure  1.  Planar  elastic  rotating  beam  problem  description. 


Figure  2.  Actuation  profiles  for  extreme  values  of  K. 
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TIME  MARCHING  METHODS 


An  implicit  Newmark  method  is  used  to  integrate  the  equations  of  motion  (Eq.  2)  in  time, 
while  a  Newton-Raphson  equilibrium  loop  within  each  time  step  handles  the  structural 
nonlinearities.  This  method  involves  a  first-order  expansion  of  the  nonlinear  internal  forces  at 
consecutive  time  steps  i  and  i+1  via  the  tangent  stiffness  matrix2: 


(P}li+i  =  (P)li  +  [K]|j  •  {Au}  {Au}  =  {u}|i+1  —  {u}|j  (3) 


Furthermore,  the  velocity  and  acceleration  terms  at  time  step  i+1  are  approximated  with  the 
Newmark  method1: 


00I.+1  = 


1 

(3  ■  At2 


■  ({Au}  —  At  •  (u}|,) 


(4) 


=  FU '  ^  -  (jp  0  ■ ■  Ml.  -  At  •  (^ - -  l)  ■  {QJI,  (5) 


where  y  and  P  are  Newmark  integration  controls,  set  to  0.5  and  0.25,  respectively.  The  three 
above  approximations  are  substituted  into  Eq.  2,  and  the  equations  of  motion  at  time  step  i+1 
then  become: 


’ [M] + FS '  [c]li«  +  [K]li) ' {Au}  =  -  tpJii 

+  [M] '  +  (Fp  “ F) ' (ii!li)  +  [c]li+1 


where  [M]  is  constant,  and  has  no  time  step  subscript.  Newton-Raphson  equilibrium  iterations 
within  each  time  step  are  undertaken  by  conditionally  updating  the  solution  {u}|i+1,  computing 
the  velocity  and  acceleration,  and  then  computing  the  residual  at  time  step  i+1: 
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{R)li+1  =  {F}|i+1  -  {P}|1+1  -  [M]  •  {u}|1+1  -  [C]|1+1  •  (u}|i+1 


(7) 


This  residual  is  then  added  to  the  external  force  vector  {F}|i+1,  {Au}  is  recomputed,  and  the 
process  repeats  until  the  residual  is  driven  to  zero. 

Having  solved  the  system  of  equations,  it  is  assumed  here  that  an  objective  function  will  take 
the  form  of  a  scalar  quantity,  integrated  in  time  : 


g=  f  G(t)  •  H(t  —  t0)  •  dt  (8) 


where  G  is  a  time-dependent  property  of  the  beam  deformation  (displacement,  stresses,  vibration 
frequency,  etc.),  and  H  is  the  Heaviside  step  function.  Thus,  it  is  only  desirable  to  optimize  the 
beam’s  performance  over  a  single  actuation  period,  once  the  motion  has  become  time -periodic, 
between  tG  and  tf.  For  reasons  that  will  be  discussed  below,  it  is  still  necessary  to  define  the 
integral  from  0  to  tf,  even  though  information  prior  to  t0  is  of  no  interest.  The  derivative  of  the 
objective  function  with  respect  to  an  unspecified  vector  of  design  variables,  {x},  is: 


dg 

3{x] 


ftf  3G  T 

l  3R 


d(u) 

3{x} 


•  H(t  —  t0)  •  dt 


(9) 


where  { u }  is  the  solution  of  the  nonlinear  system  of  equations.  Computation  of  this  derivative  is 
facilitated  by  taking  the  derivative  of  Eq.  2: 


3[M] 

3{x} 


•  {u}  +  [M]  • 


d{u}  d[C] 
3{x]  3{x} 


•  (u)  +  [C]  • 


d{u) 

3{x} 


3{u}  3{P}  _  3{F] 

L  '  3{x}  +  3{x}  3{x} 


(10) 


The  adjoint  method  is  undertaken  by  pre-multiplying  the  equations  of  motion  by  an  adjoint 
vector  {k},  and  adding  this  term  to  the  derivative  of  the  objective  function7: 
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dg 

a{x} 


/ 


-r 

“'o 


5G  T  d{u} 


a{u}  a{x} 


H(t  —  t0)  +  {A}1 


V 


/a[M]  a{u}  a[c]  \\ 


a{x} 


a{x}  a{x} 


a{u}  a{u}  a{p}  a{F} 


dt 


(11) 


Integration  by  parts  is  used  (twice)  to  remove  the  time  derivative  terms  within  the  response 
derivatives: 


dg 

a{x} 


/ 


"I 


dt 


{A}1 


^[M]  a[c]  a{p}  a{F}\ 

{uj  +  TTT  ‘  lUJ  +  -  T7  T  I  + 


a{x} 


a{x} 


a{x}  a{x}. 


aG 


\\  0{u} 


H(t  -  t0)  +  {X}T  •  [M]  -  {A}T  •  [C]  +  {Af  •  ([K]  -  [C]) 


\ 


a{u} 
a{x}  J 


(12) 


+{A}T  •  [C] 


a{u] 


a{x} 


tf 


+  (A}T  •  [Mr] 


a{u} 


a{x} 


tf  hit  fa/ii 

-  [A]  •  [M] 

o 


a{x] 


tf 


This  equation  would  suggest  the  following  boundary  conditions: 


d{u} 

a{x} 


t=o 


a{u] 


a{x} 


=  o 


t=0 


(A}lc=tf  =  {A}|t,t(  =  0 


(13) 


The  first  two  conditions  are  satisfied  due  to  the  fact  that  initial  conditions  upon  the  beam’s 
dynamics  are  prescribed  (both  position  and  velocity  are  zero  at  the  initial  condition).  If  the 
integration  bounds  were  taken  from  tG  to  tf,  rather  than  starting  at  zero,  this  would  no  longer  be 
true.  As  the  very  purpose  of  the  adjoint  method  is  to  entirely  skip  the  generally  unneeded 
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derivatives  of  the  system  response  (3{u}/3{x}),  a  direct  method  (which  does  compute  these 
terms)  would  have  to  be  marched  from  t  =  0  to  tG  in  order  to  compute  these  derivatives,  and  then 
the  adjoint  method  could  be  undertaken  from  tQ  to  tf  .  From  a  computational  cost  perspective,  it 
makes  more  sense  to  integrate  from  0  to  tf,  and  use  a  Heaviside  step  function.  The  latter  two 
relations  represent  terminal  conditions  on  the  adjoint  vector.  This  vector  can  now  be  solved  for 
with: 


T  T  •  (7  Lr 

{A}  •  [M]  —  {A}  .[C]+{A}T-([K]-[C])  =  -gj-j  •  H(t  — 10) 


(14) 


As  intended,  this  adjoint  vector  is  independent  of  the  design  variables,  and  must  be  computed 
by  integrating  backwards  in  time  with  an  implicit  Newmark  method,  starting  with  the  terminal 
conditions  and  working  backwards  to  the  initial  conditions.  Thus,  the  response  {u}  and  the 
adjoint  vector  {a}  cannot  be  computed  simultaneously.  The  system  matrices  (mass,  damping, 
and  tangential  stiffness)  at  each  time  step  then  must  either  be  stored  during  the  forward 
integration  for  future  adjoint  computations  (impractical  for  large  problems),  or  recomputed 
during  the  reverse  integration.  Having  computed  the  reduced  adjoint  vector,  the  sought-after 
design  sensitivity  is  given  by: 


dg 

5{x) 


d[M]  3[C]  5{P] 

a{x}  ‘ +  a{x}  ‘ +  a{x} 


(15) 
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MODEL  REDUCTION 


As  noted  above,  the  computational  cost  associated  with  design  optimization  via  implicit 
integration  is  large,  particularly  for  structures  with  a  large  number  of  degrees  of  freedom  (Ndof): 
each  equilibrium  iteration  within  each  time  step  requires  the  solution  to  a  set  of  equations  (Eq. 
6).  Furthermore,  the  time  step  must  be  kept  relatively  small.  This  is  not  necessarily  due  to 
stability  concerns  (the  implicit  method  is  guaranteed  stable  for  linear  systems,  though 
nonlinearities  may  lead  to  instabilities3),  but  due  to  limited  range  of  convergence  of  the  Newton- 
Raphson  algorithm:  the  initial  guess  at  time  step  i  is  the  converged  solution  at  step  i-1.  Finally, 
the  computation  of  the  adjoint  vector  (Eq.  14)  can  be  cumbersome,  as  the  terms  must  be 
reestablished  during  the  reverse  integration.  These  problems  can  all  be  attenuated,  to  some 
degree,  with  reduced  order  modeling,  wherein  the  system  of  equations  is  projected  onto  a 
reduced  basis.  This  is  done  by  approximating  the  structural  response  as  a  linear  combination  of 
basis  vectors: 


Nr 

U(t)  =  ^Pi(t)  •  4>i  {u}  =  [<b]  •  {r|}  (16) 

i=l 

where  { u }  has  been  reduced  to  a  set  of  generalized  coordinates  {q}  through  a  matrix  of  basis 
vectors,  or  modes:  [O],  The  number  of  modes  needed  (Nr)  is  typically  much  less  than  the  size  of 
the  original  system  (Ndof)-  The  reduced  basis  is  constructed  with  POD  modes,  which  requires 
computing  the  response  of  the  full  order  model  (FOM,  essentially  Eq.  2),  and  storing  the  results 
at  several  time  steps16  in  a  snapshot  matrix: 


[S]  =  [{u}li  (u}|2  -  (u}|Nsteps  ]T 


(17) 


where  [S]  is  an  Nsteps  x  NDqf  matrix.  The  following  eigenproblem  must  then  be  solved16: 


[S]T  •  [S]  •  [V]  =  [V]  •  [A] 


(18) 


The  eigenvectors  of  the  system  are  contained  in  the  matrix  [V],  and  the  eigenvalues  are  given  in 
the  diagonal  matrix  [A].  The  reduced  basis  is  then  computed  as: 


[d>]  =  [S]  •  [V] 


(19) 
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Upon  solution  of  this  eigenvalue  problem,  a  large  number  of  basis  functions  are  available,  but  a 
relatively  small  number  (Nr)  can  be  utilized  in  practice.  A  suitable  retention  number  can  be 
inferred  from  the  relative  size  of  the  eigenvalues  [A],  which  provide  the  system  energy  captured 
by  each  mode. 

The  nonlinear  system  of  equations  (Eq.  2)  can  be  projected  into  this  reduced  space  by 
introducing  the  approximate  form  of  the  beam  deformation  (Eq.  16),  and  pre-multiplying  through 
by  the  modal  matrix: 


[OjT  •  [M]  •  [cb]  •  {rj}  +  [0>]T  •  [C]  •  [O]  •  (t)}  +  [d>]T  •  {P}  =  [0]T  •  {F}  (20) 


[Mr]  •  {rj}  +  [Cr]  •  (f)}  +  {Pr}  =  {Fr} 


(21) 


where  the  subscript  (r)  implies  a  reduced  term.  The  reduced  parameters  (force  vectors  and 
stiffness,  damping,  and  mass  matrices)  can  then  be  used  with  the  above  Newmark  integration 
scheme  to  compute  {q}  at  each  time  step.  The  method  given  in  Eq.  20  requires  that  the  finite 
element  terms  be  computed  in  the  full  order  space,  assembled,  and  then  reduced  with  pre-  and 
post-multiplication.  A  superior  method,  used  here,  is  to  perform  reduction  on  an  element-by- 
element  basis,  rather  than  for  the  entire  system  .  For  example: 


[Kr] 


e=l 


[Ke]  •  [Oe] 


(22) 


where  the  global  matrix  is  found  through  a  typical  assembly  process  of  each  finite  element  (e). 
[®e]  refers  to  the  modal  matrix  terms  corresponding  to  the  degrees  of  freedom  found  in  the  eth 
finite  element,  and  is  a  subset  of  [$]. 

Using  the  reduced  order  model  alleviates  computational  cost  through  storage  and  system 
solver  requirements  (as  each  reduced  term  is  of  the  small  size  Nr  x  Nr,  albeit  full).  The  reduced 
storage  cost  also  has  significant  ramifications  for  adjoint  sensitivities,  as  discussed  below. 
Model  reduction  incurs  three  cost  penalties,  however.  First,  the  solution  of  the  POD 
eigenproblem  (Eq.  18)  can  be  costly  for  a  large  number  of  time  steps;  second,  the  assembly  of 
the  global  finite  element  matrices  (Eq.  22)  requires  more  operations  in  the  reduced  space  due  to 
the  pre-  and  post-multiplication  of  the  modal  matrix;  third  (and  typically  the  most  severe),  the 
snapshot  matrix  must  be  populated  with  data  from  the  full  order  model. 

The  following  method  is  used  to  compute  the  basis  functions.  Using  implicit  integration,  the 
full  order  model  is  integrated  from  t  =  0  to  a  specified  tpoM-  The  results  at  each  time  step  are 
stored  within  the  snapshot  matrix  (Eq.  17),  which  is  used  to  compute  the  POD  modes  (Eq.  19). 
These  modes  can  be  used  to  reduce  the  equations  of  motion,  re-integrated  from  t  =  0  to  the  last 
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time  step  of  interest  (typically  after  the  solution  becomes  time-periodic),  or  they  can  be  used  with 
the  time-periodic  methods  described  below.  The  training  period,  tpoM,  is  an  important  metric:  if 
the  value  is  too  small,  the  snapshot  matrix  will  not  adequately  capture  the  important  dynamical 
characteristics  of  the  system.  If  the  value  is  too  large,  the  cost  penalties  associated  with  running 
the  full  order  model  will  be  severe.  This  trade-off  is  clearly  demonstrated  below. 

The  sensitivities  of  the  ROM  are  found  in  a  similar  manner  to  above,  with  the  reduced 
adjoint  vector  given  as: 


T  T  •  U\Ji 

M  -  [Mr]  -  {Ar}  •[Cr]  +  {A,.}T  ([Kr]-[Cr])=-[4»]T-^  H(t-t0) 

Ml,-,,  =  ML,  =  0 


(23) 


It  should  be  noted  that,  unlike  in  the  full  order  model,  where  all  of  the  system  matrices  and 
vectors  must  be  re-established  during  the  reverse  integration,  now  all  of  the  terms  are  small 
enough  to  be  stored  in  memory  during  the  forward  integration.  As  such,  the  computational  cost 
associated  with  the  reduced  adjoint  vector  is  very  small.  It  is  also  assumed  here  that  the  reduced 
basis  (POD  modes)  is  not  a  function  of  the  design  variables,  which  is  obviously  untrue.  It  has 
been  shown  that  in  using  reduction  of  linear  structural  dynamics  to  modal  coordinates  for 
sensitivity  analysis,  the  derivative  of  the  linear  eigenvectors  with  respect  to  the  design 

o 

parameters  is  negligible  .  A  similar  assumption  with  the  POD  modes  will  be  shown  below  to  be 
acceptable.  Having  computed  the  reduced  adjoint  vector,  the  sought-after  design  sensitivity  is 
given  by: 


dg 

3{x} 


•99  + 


d[Cr] 

3{x} 


•99  + 


d{Pr) 

3{x] 


(24) 


Results  pertaining  to  ROM  accuracy,  for  both  the  system  response  and  the  adjoint 
computations,  are  given  in  Figure  3,  for  an  objective  function  g  of: 


rtf 

=  («*)' 

Jo 


H(t  —  t0)  •  dt 


(25) 


where  5tip  is  the  transverse  displacement  of  the  beam’s  tip,  measured  in  the  floating  frame  of 
Figure  1.  Minimizing  this  metric  will  essentially  decrease  the  beam  deformation  throughout  the 
actuation  cycle  (though  only  the  deformation  during  the  final,  time -periodic  motion  is  considered 
in  computing  the  objective  function).  The  normalized  tip  displacement  is  given  as  a  function  of 
time  in  the  upper  plot  of  Figure  3,  for  the  full  order  model  (integration  of  Eq.  2)  and  reduced 
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order  models  with  a  varying  number  of  retained  POD  modes  (integration  of  Eq.  21).  For  an 
actuation  frequency  of  40  rad/s  and  K  set  to  0.01,  the  response  is  relatively  nonlinear 
(approaching  50%  of  the  beam  length)  and  smooth.  The  initial  transients  decay  after  9  actuation 
cycles,  leaving  a  cyclic  sinusoidal  response  in  phase  with  the  rotation  actuation.  The  reduced 
order  models  are  computed  by  integrating  the  full  order  model  through  a  training  period  of  half 
of  an  actuation  cycle,  with  the  data  at  each  time  step  used  for  the  snapshot  matrix.  A  ROM  with 
2  modes  is  inaccurate,  diverging  from  the  correct  response  shortly  after  the  training  period,  4 
modes  has  minor  inaccuracies  but  largely  captures  the  time-periodic  response,  and  6  modes  is 
very  accurate  for  both  the  initial  transients  and  time -periodic  response. 


Figure  3.  System  response  and  adjoint  accuracy  for  various  reduced  order  models:  to  =  40 

rad/s,  K  =  0.01. 


Using  either  Eq.  14  or  Eq.  23,  the  adjoint  vector  can  be  computed  (which  is  only  effected  by 
the  objective  function  of  Eq.  25,  not  the  as-yet  unspecified  design  variables).  For  both  cases,  the 
adjoint  vector  is  analogous  to  the  system  response  due  to  an  oscillating  transverse  point  load  at 
the  tip  (equal  to  -2-  5tip)  which  only  operates  within  the  last  actuation  cycle  of  motion.  The 
member  of  the  adjoint  vector  that  corresponds  to  the  5tip  structural  degree  of  freedom  is  given  in 
the  lower  plot  of  Figure  3.  For  the  implicit  integration,  the  adjoint  has  zero  terminal  conditions 
(Eq.  13),  and  is  integrated  backwards  in  time  according  to  Eq.  14.  The  “dummy  load”8  actuates 
the  system  within  the  final  actuation  cycle,  and  is  then  removed  (as  the  objective  function  in  Eq. 
20  is  only  concerned  with  the  final  time -periodic  response).  Moving  backwards  from  9-T  to  0, 
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the  adjoints  resemble  a  damped  free  vibration  (albeit  with  a  negative  damping  matrix,  Eq.  14), 
which  decays  to  zero  as  t  approaches  0.  For  the  reduced  adjoint  vector,  a  ROM  with  2  modes  is 
expectedly  inaccurate,  though  increasing  the  number  of  modes  to  4  leads  to  instabilities  in  the 
adjoint  integration.  A  ROM  with  6  modes  is  needed  for  a  stable  and  accurate  adjoint 
computation. 

The  logarithmic  error  in  the  objective  function  and  the  norm  of  the  design  sensitivity  (Eq.  24, 
where  the  design  variables  are  the  cross-sectional  area  of  each  beam  element)  is  given  in  Figure 
4,  where  the  error  is  defined  as  the  difference  between  the  full  and  reduced  order  models.  Along 
with  the  data  of  Figure  3,  it  can  be  seen  that  more  modes  are  generally  needed  for  an  accurate 
sensitivity  analysis  than  for  the  system  response.  From  the  left  side  of  Figure  4,  there  is  an 
accuracy  trade-off  between  the  number  of  modes  needed  and  the  duration  of  the  training  period. 
If  the  training  period  is  too  short  however  (0.1-T,  for  the  example  given),  increasing  the  number 
of  modes  does  not  improve  the  accuracy.  It  can  further  be  seen  from  the  right  side  of  Figure  4 
that  the  number  of  modes  needed  for  an  accurate  ROM  is  a  weak  function  of  the  number  of 
structural  degrees  of  freedom  ,  which  should  ease  the  computational  cost  during  the  time 
integration.  The  cost  of  solving  the  full  order  system  is  expected  to  be,  in  a  best-case  scenario, 
O(Ndof)  for  each  iteration,  if  a  sparse  solver  is  employed.  The  cost  of  solving  the  reduced  order 
system  is,  conservatively,  0(Nr),  as  the  reduced  matrices  are  all  full.  Computational  cost 
savings  are  realized  if  NDqf  grows  much  faster  than  ,  as  indicated  by  Figure  4. 


Figure  4.  ROM  accuracy  as  a  function  of  the  number  of  modes  (left,  with  NDOF  fixed  at 
30),  or  the  number  of  structural  degrees  of  freedom  (right):  w  =  40  rad/s,  K  =  0,01, 
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CYCLIC-IMPLICIT  METHOD 


The  cyclic-implicit  method  uses  an  implicit  time  marching  scheme  over  a  single  actuation 
cycle  to  produce  the  relationship  between  the  system  responses  at  adjacent  time  steps,  and  then 
sets  the  system  response  at  the  final  step  to  be  identical  to  that  at  the  first  step  .  This  is  done  by 
converting  Eq.  2  into  first  order  form,  and  using  a  Crank-Nicolson  scheme  to  compute  the 
system  response  at  one  time  step  in  terms  of  the  response  at  the  previous  step1: 


(-  f[I] 

Ut  L[o] 


[0]] 
[M]  J 


[[0] 

L[o] 

\m 

'  L[o] 


-Mi 
[c]  J 
[0]1 
[M]  J 


+  (1  -  P)  • 


[[0] 

L[o] 


[A]  •  (qjli+i  +  [B]  •  {q}|j  +  P  •  {Rint}|i+i  +  (1  -  P)  -  {Rint}li 

=  P  '  (Rextlli+l  +  (1  -  P)  ’  (Rextlli 


(26) 


(27) 


where  (1  is  an  integration  parameter  (set  to  0.5  for  a  Crank-Nicolson  scheme),  and  {q},  { Rmi}. 
and  {Rext}  are  generalized  coordinates  and  forces  for  the  first  order  system.  Dividing  an 
actuation  cycle  into  Ncyc  time  steps,  and  further  assuming  that  the  solution  state  at  the  first  time 
step  and  the  final  time  step  are  identical  (for  time -periodicity),  the  following  system  of  equations 
are  obtained: 


f[A] 

[B]  [A] 

[B] 


[B] 


[A] 

[B] 
r 


- 

f  (q)  1 1 

{q}|2 

•  5 

(q}|3  > 

Jq}lNcyc  > 

R  i  n  t 

}li  +  d- 

P)  '  (Rjnt }  I  Ncyc 
P  *  (Rintlb  +  (1  -  P)  *  (Rintlll 
P  ’  (Rjnt  1 1 3  T  (1  P)  '  (Rjnt  1 1 2 


IP  '  {Rint  }lNcyc  +  (1  —  P)  •  {Rjnt }  I N cyc  -1 J 
f  P-{Rext}ll  +  (l-P)’{Rext}lNcyc  ^ 


=  < 


P-{Rext}l2  +  (l-P)'{Rext}ll 
P-{Rext}l3  +  (l-P)-{Rext}l2  2 


IP  •  {Rext  1 1 N cyc  +  (1  —  P)  •  (Rext}lNcyc -lj 


(28) 
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This  equation  can  be  re-written  as: 


[Mcyc]  •  {ucyc}  +  {Pcyc}  ~  {pcyc} 


djPcyc] 

a{ucyc} 


(29) 


Where  { Ucyc } ,  { Pcyc }  -  and  { Fcyc }  are  displacements/velocities  and  forces  grouped  by  common 
temporal  degree  of  freedom,  for  each  time  step  in  the  cycle.  Solution  of  this  nonlinear  system  is 
facilitated  with  a  nonlinear  Jacobian  matrix  and  a  residual  vector: 


[jcyc]  [^cyc]  T  [KCyC  ] 


(30) 


{Rcyc  }  -  [Mcyc  ]  '  {UCyc  }  +  {Pcyc  }  _  h  ‘  {pcyc  } 


(31) 


The  system  can  be  solved  with  a  Newton-Raphson  equilibrium  loop,  gradually  increasing  p  (load 
scaling  factor)  from  0  to  1  to  ensure  convergence  of  the  solution.  The  solution  guess  at  the  nlh 
iteration  is: 


{Ucyc  1  ={u 

eye  }  [j  eye  p.{R 

eye  } 


(32) 


This  process  is  repeated  until  the  residual  vector  vanishes,  at  which  point  convergence  has  been 
obtained. 

For  the  cyclic-implicit  method,  the  time-integrated  objective  function  is  defined  with  a 
numerical  integration: 


geye  =  {^}T  •  {Gcyc}  {Gcyc}  =  {G|i  G|2  •••  G|Ncyc]T  (33) 


where  {ro}  is  a  vector  of  integration  weights.  Adjoints  are  computed  by  first  differentiating  the 
objective  function: 


dgeye  _  dgeye  T  d{UCyc}  _  T  .  3(Gcyc }  g{Ucyc} 

a{x]  a{ucyc}  a{x}  “  a{ucyc}  a{x} 


15 


The  derivative  of  the  system  of  equations  (Eq.  29)  is: 


^.^cyc.  (-  1  I  Tt  1  5^Ucyc }  I  ^{^cyc}  _ 

d{x )  ‘  lUcyc)  +  LJcycJ  •  a{x)  +  a{x}  d{x} 


(35) 


As  above,  this  equation  is  pre-multiplied  by  an  adjoint  vector,  and  added  to  the  derivative  of  the 
objective  function: 


d&cyc 

3{x} 


d{GCyc}  d{UCyc}  -iT 

a{ucyc}'  d{x)  +[Acyc] 

/a[ Mcyc]  r  I  r  1 
l  d{x}  lucyc  j  '  Ltcyc  J 


g{Ucyc} 

3{x} 


^{Pcyc  ) 
5{x} 


^{Fcyc}\ 

a{x}  J 


(36) 


This  system  of  equations  for  the  solution  of  the  adjoint  vector  then  becomes  simply: 


{\yc  }T  •  [jcyc  ]  =  -{w}T  • 


d{Gcyc } 

3{Ucyc) 


(37) 


Finally,  the  sought-after  design  sensitivities  are  given  by: 


^{Fcyc}\ 

a{x}  J 


(38) 


The  cyclic-implicit  method  can  be  used  in  either  the  full  or  the  reduced  order  space,  by 
converting  either  Eq.  2  or  Eq.  21  into  the  first  order  form  of  Eq.  27.  In  the  event  that  the  ROM  is 
to  be  used,  time  marching  of  the  full  order  model  must  be  run  through  the  training  period  to 
compute  the  POD  modes,  as  discussed  above. 
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SPECTRAL  ELEMENT  METHOD 


The  spectral  element  method  is  a  monolithic  time  approximation,  best  described  by  first 
considering  a  nonlinear  system  with  a  single  degree  of  freedom,  u.  The  spectral  element  method 
discretizes  the  time  domain  into  a  finite  number  of  elements.  The  width  of  the  jth  spectral 
element  (in  seconds)  is  hJ.  Each  spectral  element  is  transformed  onto  the  interval  (  E  [-1,  1],  and 
the  solution  of  the  dependent  variable  within  the  element  is  approximated  by  an  mth-order 
Lagrange  polynomial"  : 


fl'  =  ^  ui&)  •  <© 

k=0 


(39) 


where  u'  is  the  approximate  solution,  uJ((k)  is  the  solution  (unknown  nodal  degrees  of  freedom  in 
time)  sampled  at  m  locations  through  the  element  (<^  are  the  m-1  zeros  of  the  Lobatto 
polynomials,  as  well  as  -1  and  1),  and  i|/kj  is  the  kth  Lagrange  polynomial  within  spectral  element 
j.  A  non-dimensional,  second-order,  single-degree-of-freedom  system  is  defined  within  the  jth 
spectral  element  by: 


4  •  M 

(W 


d2ui  2-C 
d^“  +  _hT 


duj 


•  —  +  pJ(0  =  f,(0 


(40) 


where,  as  above,  M  is  the  mass,  C  the  damping,  p  is  the  internal  force  (a  nonlinear  function  of  u), 
and  f  is  the  external  force.  The  approximate  solution  given  above  can  be  substituted  into  this 
equation,  and  the  residual  is  minimized  with  the  Bubnov -Galerkin  method4.  The  kth-order 
Lagrange  polynomial  is  used  as  a  trial  function: 


d2uj  duj  hj  /  .  \\ 

+  c  •  +  y  (p)  (0  -  f)  (0)  )  *  d?  =  0 


d?: 


(41) 


Integration  by  parts  (twice)  transfers  the  time  derivatives  from  the  approximate  solution  to  the 
trial  function: 


2  •  M 
hi 


+  C  ■  •  O' 


1 

+ 

-1 


(42) 
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•  d^  =  0 


2  •  M 
hi 


dX  r  £^L\ 

d^2  d^  J 


hj  ,  . 

+  y.(p,-f)) 


The  approximate  solution  can  be  expanded  in  terms  of  the  nodal  degrees  of  freedom  at  the  m 
Lobatto  points  (Eq.  39),  m  trial  functions  are  considered,  and  the  integral  terms  are  approximated 
with  a  quadrature  rule: 


1  m 

[  S(0-d?  =  ^S(?k)-a)k 

-1  k=0 


(43) 


where  S  is  any  function  of  ^  and  is  the  Gaussian  quadrature  weight  at  Lobatto  point  k.  The 
equation  of  motion  is  now  given  by  the  following  system  of  equations  within  the  jth  spectral 
element: 


M  •  K]  •  {uU}  +  C  ■  [v|]  •  {u'se}  +  [iL]  ■  {p'se}  =  [iL]  ■  {fie}  (44) 


Where  [Tm]  and  [*TC]  are  matrices  of  2nd  and  1st  derivatives  of  the  Lagrange  polynomials 
evaluated  at  the  Lobatto  points,  and  [Im]  is  a  diagonal  matrix  containing  the  Lobatto  weights. 
Explicit  formulations  are  provided  in  Ref.  22.  The  vectors  in  Eq.  44  are  a  collocation  of  terms  at 
each  node  within  the  jth  spectral  element. 

The  solution  over  the  entire  domain  can  be  obtained  with  a  typical  assembly  process  of  the 
spectral  elements,  adding  contributions  to  shared  nodes  (namely  the  solution  at  ^  =  +  1  within 
each  element)  into  the  same  location  within  the  global  matrices.  Lurthermore,  time -periodic 
solutions  can  be  obtained  by  assuming  that  <^0  of  the  first  element  and  (Jm  of  the  last  element  are 
the  same  node.  Contributions  within  the  global  matrices  pertaining  to  this  last  node  are  added  to 
the  contributions  pertaining  to  the  first  node4.  The  size  of  the  assembled  system  is  thus  reduced 
from  (m  x  Nei)  +  1  to  m  x  Nei,  where  Nei  is  the  number  of  spectral  elements,  and  m  is  the  order  of 
the  approximation  made  within  each  element.  Lor  cases  with  more  than  one  structural  degree  of 
freedom,  the  final  system  of  equations  can  be  computed  as: 

([M]  ®  [Vm]  +  [C]  ®  [%])  •  {Use}  +  [I]  ®  [I J  •  {Pse} 

=  [I]  0  [IJ  ■  (Fse)  [Kse]  =  (45) 


where  (x)  denotes  a  tensor  product,  and  [I]  is  the  identity  matrix  of  size  NdofX  NDof,  where  NDof 
is  the  size  of  the  finite  element  system.  The  vectors  { XJse } ,  { Pse } ,  and  { Lse } ,  are  a  collocation  of 
the  deformation,  internal,  and  external  forces,  respectively,  at  each  nodal  time,  grouped  by 
common  structural  degrees  of  freedom.  Lor  example: 
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(Use)  =  < 
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UnD0F  ltl 
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>  ...  < 
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> 

V 
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UNDof  L 

^  tm  xNei  y 

> 

(46) 


Similar  quantities  are  given  for  the  nonlinear  internal  force  and  the  external  force.  As  above,  the 
solution  of  this  system  is  facilitated  with  a  nonlinear  Jacobian  matrix  and  a  residual  vector: 


[Jse ]  =  [M]  ®  [Vm]  +  [C]  ®  [¥c]  +  [I]  ®  [Ij  •  [Kse] 


(47) 


(Rse)  =  ([M]  ®  [Vm]  +  [C]  ®  [%])  •  {Use}  +  [I]  ®  [I  J  •  ({Pse}  -  p  •  {Fse}) 


(48) 


The  system  can  be  solved  with  a  Newton-Raphson  equilibrium  loop,  as  in  Eq.  32. 

A  sensitivity  analysis  of  the  spectral  element  method  proceeds  in  a  very  similar  manner  to  the 
above  cyclic-implicit  method.  The  objective  function  is  defined  with  a  quadrature  rule: 


gse  =  !<'>1T  •  f !’•'>)  =  diagdXJ)  =  {Git,  Gltz 


Gl.„»d}T  («) 


This  system  of  equations  for  solution  of  the  adjoint  vector  is: 


{^se}T  •  [Jse]  =  -{w}T  • 


d(Gse} 

3{Use} 


(50) 


The  design  sensitivities  are  given  by: 


=  Re}T 


fd([M]  ®  [%,])  riT  ,  ,  a([C]  ®  [Vc]) 

V  3{x}  ‘1Usej+  5{x} 


•{Use}  +  [I]0[IJ 


djPse)  d{¥se}\\ 

a{x}  a{x}  ) J 


(51) 
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As  above,  the  spectral  element  techniques  given  in  this  section  can  be  used  for  either  the  full  or 
the  reduced  order  modeling,  by  substituting  the  appropriate  reduced  terms  in  the  various 
equations. 

A  comparison  between  the  spectral  element  method  and  the  cyclic -implicit  method  is  given 
in  Figure  5;  the  logarithmic  error  in  the  objective  function  and  the  norm  of  the  gradient  vector  is 
plotted  as  a  function  of  the  number  of  time  steps  per  actuation  cycle.  The  error  is  given  as  the 
difference  between  the  current  solution  and  the  solution  computed  with  600  temporal  degrees  of 
freedom.  For  the  spectral  element  results,  the  order  of  the  Lagrange  basis  functions  (Eq.  39)  is 
fixed  at  m  =  5;  the  temporal  degrees  of  freedom  are  expanded  by  adding  more  elements  (h- 
refinement).  As  would  be  expected,  the  higher-ordered  spectral  elements  are  more  accurate 
(lower  error)  than  the  cyclic-implicit  method. 


temporal  degrees  of  freedom 


Figure  5.  Error  metrics,  with  the  polynomial  order  of  the  spectral  elements  (m)  fixed  at  5: 

NDOF  =  30,  o>  =  15  rad/s,  K  =  0.99. 

Both  time -periodic  methods  have  the  option  of  grouping  the  solution  vector  ( { Ucyc }  or  {Use}) 
by  common  structural  degrees  of  freedom  (as  formulated  for  the  spectral  element  method,  Eq. 
46)  or  by  common  temporal  degrees  of  freedom  (as  formulated  for  the  cyclic-implicit  method, 
Eq.  28).  The  two  choices  give  drastically  different  sparsity  patterns  in  the  Jacobian  matrices,  as 
seen  in  Figure  6.  For  the  spectral  element  method,  grouping  by  common  structural  degree  of 
freedom  is  preferable,  as  the  matrix  is  closely  banded  towards  the  diagonal.  For  the  cyclic- 
implicit  method,  the  matrices  are  twice  as  large,  due  to  the  first  order  implementation.  The 
method  can  of  course  be  formulated  as  a  second  order  system,  but  this  would  require  a  three  time 
level  scheme,  rather  than  the  two  time  levels  seen  in  Eq.  27,  complicating  the  process.  Grouping 
by  common  temporal  degree  of  freedom  is  slightly  superior  for  the  cyclic-implicit  method, 
though  neither  skyline  is  very  efficient  from  a  sparse-solver  standpoint1. 
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Figure  6.  Jacobian  sparsity  patterns:  NDOF  =  30,  to  =  15  rad/s,  K  =  0.99. 


Neither  Figure  5  nor  Figure  6  provides  an  impetus  for  the  use  of  the  cyclic-implicit  method 
over  the  spectral  element  method:  it  isn’t  as  accurate,  the  Jacobians  are  twice  as  large,  and  the 
computational  effort  required  for  a  sparse  solver  is  much  higher.  As  will  be  seen  below, 
however,  the  cyclic-implicit  method  tends  to  drive  the  residual  of  the  nonlinear  system  of 
equations  (Eq.  31)  to  zero  at  a  faster  rate,  without  relying  as  heavily  upon  the  load  scaling  factor. 
This  becomes  particularly  important  for  strong  nonlinearities,  where  the  spectral  element  system 
of  equations  (Eq.  45)  may  require  the  load  scaling  factor  to  be  slowly  increased  from  0  to  1  to 
prevent  the  system  from  diverging. 
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OPTIMIZATION  RESULTS 


Using  the  6  techniques  derived  above  (implicit  Newmark  time  marching,  spectral  elements, 
and  the  cyclic-implicit  method,  each  with  or  without  POD-based  model  reduction),  it  is  desired 
to  minimize  the  cycle-integrated  tip  deflection  (Eq.  25),  using  the  cross-sectional  area  of  each 
beam  element  as  the  deign  variables  { x } .  No  constraints  are  considered,  other  than  side 
constraints.  Each  element’s  cross-sectional  area  is  required  to  be  between  2  cm”  and  6  cm".  The 
initial  design  is  a  uniform  beam  with  an  area  of  4  cm2.  Gradient-based  optimization  is  run  with  a 
simple  steepest-descent  algorithm. 

For  smooth  kinematics,  decreasing  the  design  parameters  { x }  will  decrease  the  stiffness  of 
the  beam,  increasing  the  deformation.  Increasing  this  variable  will  increase  the  weight  and  thus 
the  inertial  loads,  also  increasing  the  deformation:  an  optimal  design  will  utilize  a  compromise 
between  the  two  design  philosophies.  This  can  be  seen  in  Figure  7,  where  (for  sinusoidal 
kinematics  with  K  =  0.01)  the  optimal  wing  structure  is  a  tapered  beam.  The  beam  is  very  stiff  at 
the  root  (where  the  bending  moments  are  largest),  and  very  thin  at  the  tip,  where  the  inertial 
loads  are  largest.  The  side  constraints  upon  the  cross-sectional  area  are  enforced  at  both  the  root 
and  the  tip.  Linear  and  nonlinear  deformations  require  the  same  beam  geometry,  as  the  designs 
given  at  15  and  40  rad/s  are  identical.  These  findings  are  not  the  case  for  pseudo-triangular 
actuation  motions  (K  =  0.99).  Deviations  from  the  baseline  are  relatively  small  (potentially 
indicative  of  a  local  minimum  within  the  design  space),  and  do  not  approach  the  side  constraints. 
Contrasting  design  philosophies  are  required  for  15  and  40  rad/s,  with  the  former  advocating 
slightly  increasing  the  beam  size  towards  the  wing  tip. 


Figure  7.  Optimal  cross-sectional  area  distributions  for  different  actuation  frequencies  and 

kinematics,  with  ten  finite  elements. 


The  deformation  of  the  beam’s  tip  (5tip)  is  given  in  Figure  8,  for  a  sinusoidal  actuation  (K  = 
0.01)  at  15  rad/s.  Both  the  transient  solution  (computed  through  implicit  integration  via  the 
Newmark  method)  and  the  periodic  solution  (computed  with  spectral  elements  and  the  cyclic- 
implicit  method,  and  compared  to  the  last  cycle  of  the  Newmark  method)  are  given,  computed 
with  the  full  order  model.  For  the  implicit  integration,  6  cycles  are  needed  for  the  deformation  to 
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reach  a  time -periodic  response,  and  the  solution  is  computed  at  100  time  steps  per  cycle.  Four 
spectral  elements  (Nei)  are  utilized  with  a  5  th  order  polynomial  (m)  approximation  within  each 
element,  while  20  time  steps  (Ncyc)  are  used  for  the  cyclic-implicit  method.  It  can  be  seen  that 
the  optimal  design  given  in  Figure  8  decreases  the  tip  displacement  from  4%  to  1%  of  the  beam’s 
length.  The  simple  deformation  patterns  (the  natural  frequency  of  the  beam  itself  is  evident  in 
the  initial  transients,  but  damps  out:  the  remaining  response  is  a  sinusoidal  trend  at  the  actuation 
frequency)  are  well-captured  by  both  time-periodic  methods  for  the  baseline  and  the  optimal 
designs.  Only  data  from  the  full  order  models  are  given  in  Figure  8:  reduced  order  modeling 
data  should  be  nearly  identical,  provided  the  number  of  modes  and  the  length  of  the  training 
period  is  sufficient  (Figure  4). 


transient  solution 


0  2  4  6 

t/T 

Figure  8.  Transient  (left)  and  periodic  (right)  beam  tip  deflection  for  the  baseline  and 
optimal  designs,  computed  with  the  full  order  model:  to  =  15  rad/s,  K  =  0.01. 


periodic  solution 


The  computational  cost  associated  with  the  optimization  of  Figure  8  is  given  in  Figure  9,  for 
both  10  and  200  beam  elements.  Such  a  fine  discretization  is  not  needed  for  this  problem,  as 
both  cases  give  identical  results,  but  the  study  will  give  an  idea  of  computational  speed-up  for 
future  cases  with  high  degrees  of  freedom.  The  computational  cost  of  each  design  iteration  is 
given  for  all  six  methods,  and  includes  computer  time  needed  for  both  the  system  response 
computation  and  the  adjoint  sensitivities,  and  in  the  case  of  the  3  reduced  models,  the  implicit 
integration  of  the  full  order  model  through  the  training  period  for  POD  mode  computations.  For 
this  case,  5  modes  were  used  for  10  beam  elements,  7  modes  for  200  elements,  both  with  a 
training  period  of  0.1-T.  All  six  methods  provide  nearly-identical  responses  at  each  iteration  and 
converge  to  the  same  design,  though  computational  effort  varies. 

The  implicit  integration  of  the  full  order  model  has  the  highest  cost,  while  model  reduction  of 
the  Newmark  method  provides  a  speed-up  by  a  factor  of  2.  Roughly  half  of  this  speed-up  is  due 
to  the  reduced  cost  of  the  adjoint  vector  (Eq.  23,  where  the  reduced  terms  are  stored  in  memory 
during  the  forward  integration);  the  remainder  is  due  to  reduced  solver  costs.  The  relative  speed- 
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up  is  expected  to  increase  as  the  size  of  the  full  order  model  (Ndof)  increases  ,  as  well  as  with 
higher  geometric  complexities  (with  less-favorable  sparse  skylines  than  the  current  case). 


_4  10  elements:  NL_„  =  30  200  elements:  bL__  =  600 

1^  DOF  DOF 


Newmark,  FOM 
Newmark,  ROM 
SE,  FOM 
SE,  ROM 
eye.,  FOM 
eye.,  ROM 


computer  time  [si 


computer  time  [si 


Figure  9.  Optimization  cost  of  full  and  reduced  order  models:  to  =  15  rad/s,  K  =  0.01. 


For  this  particular  case  (smooth  sinusoidal  kinematics  with  relatively  linear  structural 
dynamics),  the  time-periodic  methods  are  very  fast.  Convergence  of  the  system  is 
straightforward  (the  nonlinear  systems  can  be  solved  directly  with  p  =  1),  and  relatively  few 
temporal  degrees  of  freedom  are  needed.  The  time -periodic  methods  provide  a  speed-up  factor 
of  20,  with  no  significant  differences  between  the  spectral  element  method  and  the  cyclic- 
implicit  method  for  10  beam  elements.  The  cost  alleviation  through  model  reduction  (namely, 
smaller  Jacobians)  is  offset  by  the  cost  of  the  training  period  for  the  POD  modes;  model 
reduction  of  the  time-periodic  methods  has  few  benefits  for  this  case.  Contrastingly,  for  the 
larger  problem  (200  elements)  the  full  order  cyclic-implicit  method  becomes  very  slow  (due  to 
the  poor  skyline  of  the  Jacobian,  Figure  6),  and  model  reduction  can  be  used  to  drastically 
decrease  the  design  cost.  Increasing  the  discretization  from  10  to  200  elements  provides  no 
significant  differences  in  the  relative  cost  associated  with  the  remaining  5  methods. 

Similar  results  are  given  in  Figure  10  and  Figure  11,  for  sinusoidal  actuation  at  40  rad/s, 
where  the  system  response  of  the  full  order  model  is  reproduced  from  Figure  3.  The  higher 
actuation  frequency  provides  a  fairly  nonlinear  response:  the  amplitude  of  the  baseline  design’s 
deformation  is  increased  from  4%  to  40%  of  the  beam’s  length,  though  optimization  is  able  to 
drop  the  displacement  to  10%.  Ten  actuation  cycles  are  required  for  the  implicit  integration  to 
settle  within  a  periodic  response,  with  200  steps  per  cycle.  As  before,  the  harmonic  content  of 
the  time -periodic  response  is  low:  only  4  spectral  elements  are  needed  with  a  5th  order 
approximation  within  each  element,  while  20  time  steps  are  used  for  the  cyclic-implicit  method. 
Ten  POD  modes  are  used  to  build  the  reduced  models,  with  a  training  period  of  0.2- T.  The 
relative  cost  of  each  method  is  similar  to  the  previous  case,  despite  the  stronger  nonlinearities, 
with  the  cyclic-implicit  method  showing  a  very  high  cost  with  finer  discretization. 
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transient  solution  periodic  solution 


Figure  10.  Transient  (left)  and  periodic  (right)  beam  tip  deflection  for  the  baseline  and 
optimal  designs,  computed  with  the  full  order  model:  to  =  40  rad/s,  K  =  0.01. 
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Figure  11.  Optimization  cost  of  full  and  reduced  order  models:  to  =  40  rad/s,  K  =  0.01. 

Results  from  the  pseudo-triangular  actuation  kinematics  (K  =  0.99,  Figure  2)  are  given  in 
Figure  12  and  Figure  13,  for  a  frequency  of  15  rad/s.  As  expected,  a  very  high  harmonic  content 
can  be  seen  in  the  structural  dynamics  within  each  cycle.  The  inertial  forces  are  essentially  zero 
for  the  linear  portions  of  the  kinematics  (seen  in  Figure  2),  while  a  sudden  nonzero  acceleration 
is  introduced  at  the  top  and  bottom  of  each  stroke.  For  the  implicit  integration,  6  cycles  are 
needed  to  reach  periodicity,  with  300  steps  per  cycle.  For  the  spectral  element  method,  25 
elements  are  needed,  with  a  5th  order  approximation  in  each.  Presumably,  elements  could  be 
grouped  towards  the  stroke  reversal  portions  of  the  actuation  for  higher  accuracy4,  though  this  is 
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not  done  here.  For  the  cyclic-implicit  method,  200  time  steps  are  needed  during  the  actuation, 
nearly  twice  the  temporal  degrees  of  freedom  required  for  the  spectral  element  method,  as  seen 
in  Figure  5.  As  above,  the  time -periodic  and  the  Newmark  implicit  integration  are  in  close 
agreement  for  both  the  baseline  and  the  optimal  designs,  both  of  which  show  a  bi-harmonic 
response,  with  the  vibration  frequency  roughly  4  times  faster  than  the  actuation  frequency. 
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Figure  12.  Transient  (left)  and  periodic  (right)  beam  tip  deflection  for  the  baseline  and 
optimal  designs,  computed  with  the  full  order  model:  to  =  15  rad/s,  K  =  0.99. 
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Figure  13.  Optimization  cost  of  full  and  reduced  order  models:  0)  =  15  rad/s,  K  =  0.99. 

Computational  costs  (Figure  13)  show  a  significant  redistribution  of  the  run  times  of  the  6 
methods  (where  each  ROM  uses  6  and  12  POD  modes  for  the  coarse  and  fine  discretization, 
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respectively):  specifically,  the  speed-up  of  the  time-periodic  methods  becomes  less  pronounced, 
as  the  Jacobians  become  very  large  with  high  temporal  degrees  of  freedom,  and  system  solver 
costs  escalate.  The  full  order  cyclic-implicit  method  becomes  very  uncompetitive,  with  a 
computational  cost  3  times  higher  than  the  full  order  Newmark  integration.  Model  reduction  can 
be  used  to  alleviate  this  cost,  of  course,  with  the  two  reduced  time -periodic  methods  providing 
the  lowest  design  cost  for  this  case.  The  speed-up  of  the  reduced  spectral  element  method  over 
the  full  method  is  less  drastic  than  seen  with  the  cyclic-implicit  method,  as  the  full  Jacobians  are 
well-banded  (Figure  6). 

Finally,  results  are  given  in  Figure  14  and  Figure  15  for  an  actuation  frequency  of  40  rad/s 
and  a  K  of  0.99.  The  beam  vibrates  at  roughly  twice  the  actuation  frequency,  with  a  strongly 
nonlinear  response  (40%  of  the  beam  length,  with  the  optimal  design  slightly  less,  at  30%) :  this 
presents  a  very  difficult  case  for  the  time-periodic  solvers  to  compute.  For  the  implicit 
integration,  10  cycles  are  needed  to  reach  periodicity,  with  200  steps  per  cycle.  For  the  spectral 
element  method,  25  elements  are  needed  with  a  5th  order  approximation  in  each,  while  125  time 
steps  are  needed  for  the  cyclic -implicit  method.  The  efficacy  of  the  spectral  element  method 
(either  full  or  reduced)  suffers  for  this  case,  as  a  fairly  low  load  scaling  factor  (gradually 
increased  to  1)  is  needed  to  prevent  the  nonlinear  system  of  equations  (Eq.  45)  from  diverging. 
As  noted  above,  the  cyclic-implicit  method  shows  superior  convergence  characteristics  within 
the  Newton-Raphson  loop,  and  has  the  lowest  computational  cost  of  any  of  the  choices  in  Figure 
15  (as  long  as  model  reduction  is  used).  It  can  also  be  seen  than  the  reduced  Newmark 
integration  has  relatively  low  design  costs  for  this  case:  time -periodic  methods  become  less 
favorable  for  structural  dynamics  cases  with  high  frequency  content  and  strong  nonlinearities. 


transient  solution  periodic  solution 


Figure  14.  Transient  (left)  and  periodic  (right)  beam  tip  deflection  for  the  baseline  and 
optimal  designs,  computed  with  the  full  order  model:  to  =  40  rad/s,  K  =  0.99. 
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Figure  15.  Optimization  cost  of  full  and  reduced  order  models:  to  =  40  rad/s,  K  =  0.99. 
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CONCLUSIONS 


This  work  has  detailed  efforts  intended  to  reduce  the  computational  cost  associated  with  the 
gradient-based  design  of  nonlinear  rotary  beam  vibrations.  The  equations  of  motions,  formulated 
with  a  nonlinear  corotational  method  and  standard  multibody  dynamics  techniques,  are  solved 
with  six  methods:  implicit  time  marching,  a  cyclic-implicit  method,  and  a  spectral  element 
method,  each  with  or  without  POD-based  model  reduction.  Upon  computation  of  the  system 
response,  gradients  of  the  transient  response  with  respect  to  a  large  vector  of  structural  design 
parameters  are  computed  with  the  adjoint  method,  to  obtain  a  computational  cost  nearly 
independent  of  the  number  of  design  variables.  This  information  is  then  used  for  gradient-based 
optimization.  The  following  conclusions  can  be  drawn: 

1.  For  implicit  integration,  the  low  storage  costs  associated  with  model  reduction  allow  for 
extremely  inexpensive  reverse  integration  of  the  adjoint  vector,  though  more  modes  are 
typically  needed  for  an  accurate  adjoint  computation  than  for  the  forward  integration  of  the 
system  response.  Stability  can  be  a  concern  during  the  adjoint  computation  if  the  number  of 
modes  is  insufficient. 

2.  The  addition  of  more  modes  to  an  implicitly  integrated  reduced  order  model  allows  for  a 
shorter  training  period,  to  an  extent.  The  number  of  modes  required  for  a  given  level  of 
accuracy  is  a  weak  function  of  the  number  of  structural  degrees  of  freedom. 

3.  The  time-periodic  spectral  element  solver  is  typically  more  accurate  than  the  cyclic-implicit 
method  for  a  given  number  of  temporal  degrees  of  freedom,  and  also  has  a  much  more 
benevolent  skyline  for  solving  sparse  systems  of  equations.  The  latter  method,  however,  has 
better  convergence  characteristics  within  a  nonlinear  Newton-Raphson  loop. 

4.  Due  to  unfavorable  sparsity  patterns,  the  computational  cost  of  the  cyclic -implicit  method 
becomes  very  large  for  problems  requiring  high  structural  or  temporal  degrees  of  freedom; 
POD-based  model  reduction  can  be  used  to  attenuate  this  cost  substantially. 

5.  Implicit  integration  with  model  reduction  consistently  provides  a  speed-up  factor  between  2 
and  2.5  over  full  order  time  marching  for  the  problem  considered  here.  This  number  is 
expected  to  improve  for  structures  with  greater  geometric  complexity  requiring  finer  meshes. 

6.  For  problems  with  a  low  harmonic  content  (smooth  sinusoidal  actuation),  either  time- 
periodic  method,  with  or  without  model  reduction  (with  the  exception  given  in  point  4),  will 
substantially  improve  the  design  cost  over  time-marching  methods. 

7.  For  relatively  linear  problems  with  a  high  harmonic  content  (i.e.,  requiring  many  temporal 
degrees  of  freedom),  model  reduction  must  be  used  to  limit  the  computational  cost  of  the 
time-periodic  methods.  For  stronger  nonlinearities,  the  spectral  element  method  suffers  due 
to  convergence  problems  within  the  Newton-Raphson  loop:  the  reduced  cyclic-implicit 
method  is  shown  to  be  superior. 
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